!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   FEMilaro -- Finite Element Method toolkit
!!
!! FEMilaro is free software; you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 3 of the License, or
!! (at your option) any later version.
!!
!! FEMilaro is distributed in the hope that it will be useful, but
!! WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
!! General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with FEMilaro; If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>

!> Test the \c mod_eigen module.
program eigen_test

 use mod_messages, only: &
   mod_messages_constructor, &
   mod_messages_destructor, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_constructor, &
   mod_kinds_destructor, &
   wp

 use mod_eigen, only: &
   mod_eigen_constructor, &
   mod_eigen_destructor,  &
   eigen_nonsymm

 implicit none

 integer :: ne
 real(wp) :: er(8), ei(8), vr(8,8), vi(8,8)
 ! define a general square matrix a
 real(wp) :: a(8,8), m(8,8)

 a = reshape((/            &
  -9.0_wp,-3.0_wp,  6.0_wp,-3.0_wp,-6.0_wp, 9.0_wp,-10.0_wp, 6.0_wp, &
   5.0_wp, 7.0_wp, -1.0_wp,10.0_wp,-3.0_wp, 6.0_wp, -6.0_wp, 7.0_wp, &
  -2.0_wp, 6.0_wp, -4.0_wp, 3.0_wp,-2.0_wp,-4.0_wp,  3.0_wp,-1.0_wp, &
  -8.0_wp, 0.0_wp, -7.0_wp, 2.0_wp,-7.0_wp, 4.0_wp,  3.0_wp, 9.0_wp, &
  -6.0_wp, 1.0_wp, -1.0_wp, 7.0_wp, 4.0_wp,-8.0_wp, -9.0_wp,-2.0_wp, &
   0.0_wp, 3.0_wp,-10.0_wp,-6.0_wp, 3.0_wp,-8.0_wp, -4.0_wp,-9.0_wp, &
   9.0_wp, 6.0_wp, -6.0_wp, 8.0_wp, 5.0_wp,-3.0_wp,  2.0_wp, 1.0_wp, &
  -6.0_wp,-4.0_wp, -9.0_wp,-0.0_wp, 8.0_wp, 4.0_wp,  9.0_wp,-2.0_wp/),&
             (/8,8/),order=(/2,1/))

 call mod_messages_constructor()
 call mod_kinds_constructor()
 call mod_eigen_constructor()

 m = a
 write(*,'(8e14.6)') transpose(m)
! call eigen_nonsymm(m,ne,er,ei,vr,vi)
 write(*,*) 'Schur form'
 write(*,'(8e14.6)') transpose(m)
 write(*,*) 'Eigenvalues'
 write(*,'(8e14.6)') er
 write(*,'(8e14.6)') ei

 ! Now we make an experiment with a matrix which is not reduced
 m = a
 m(4:,1:3) = 0.0_wp
 write(*,*) ''
 write(*,'(8e14.6)') transpose(m)
! call eigen_nonsymm(m,ne,er,ei,vr,vi)
 write(*,*) 'Schur form'
 write(*,'(8e14.6)') transpose(m)
 write(*,*) 'Eigenvalues'
 write(*,'(8e14.6)') er
 write(*,'(8e14.6)') ei

 ! In this case, then, some eigenvalues are already evident in the
 ! upper left corner
 m = a
 m(2:,1) = 0.0_wp
 m(4:,2:3) = 0.0_wp
 write(*,*) ''
 write(*,'(8e14.6)') transpose(m)
! call eigen_nonsymm(m,ne,er,ei,vr,vi)
 write(*,*) 'Schur form'
 write(*,'(8e14.6)') transpose(m)
 write(*,*) 'Eigenvalues'
 write(*,'(8e14.6)') er
 write(*,'(8e14.6)') ei

 ! Now a symmetric matrix: the Schur form must be diagonal
 m = 0.5_wp*(a+transpose(a))
 write(*,*) ''
 write(*,'(8e14.6)') transpose(m)
 call eigen_nonsymm(m,ne,er,ei,vr,vi)
 write(*,*) 'Schur form'
 write(*,'(8e14.6)') transpose(m)
 write(*,*) 'Eigenvalues'
 write(*,'(8e14.6)') er
 write(*,'(8e14.6)') ei
 write(*,*) 'Eigenvectors'
 write(*,'(8e14.6)') transpose(vr)
 write(*,*)
 write(*,'(8e14.6)') transpose(vi)
 write(*,*) 'Residuals'
 m = 0.5_wp*(a+transpose(a))
 do ne=1,size(vr,2)
   write(*,*) hypot( &
     norm2(matmul(m,vr(:,ne))-(er(ne)*vr(:,ne)-ei(ne)*vi(:,ne))) , &
     norm2(matmul(m,vi(:,ne))-(er(ne)*vi(:,ne)+ei(ne)*vr(:,ne))) )
 enddo

 call mod_eigen_destructor()
 call mod_kinds_destructor()
 call mod_messages_destructor()

end program eigen_test

